function V = slpoolvar(X, w, vmean)
%SLPOOLVAR Computes the pooled variance
%
% [ Syntax ]
%   - V = slpoolvar(X, w)
%   - V = slpoolvar(X, w, vmean)
%
% [ Arguments ]
%   - X:            The sample matrix (d x n)
%   - w:            The matrix of weights (K x n)
%   - vmean:        The matrix of mean vectors (d x K)
%   - V:            The resultant vector of pooled variances (d x 1)
%
% [ Description ]
%   - V = slpoolvar(X, w) computes the pooled variances defined as:
%      $ v(k) = (1/sum_{ij} w)\sum_i\sum_j w_{ij} (x_{ki} - \mu_{kj})^2 $
%
%     Suppose there are n samples in d-dimensional space, and K models, 
%     then X should be a d x n matrix, and w should be a K x n matrix
%     of weights. 
%
%   - V = slpoolvar(X, w, vmean) computes the pooled variance vector
%     with the mean vectors given. 
%
%     Suppose, there are K models, then vmean should be a matrix of size
%     d x K, with each column giving the mean vector of a particular 
%     model.
%
% [ History ]
%   - Created by Dahua Lin, on Dec 24, 2007
%

%% parse and verify input

error(nargchk(2, 3, nargin));

assert(isnumeric(X) && ndims(X) == 2, ...
    'sltoolbox:slpoolcov:invalidarg', ...
    'X should be a numeric matrix.');

[d, n] = size(X);

assert(isnumeric(w) && ndims(w) == 2 && size(w,2) == n, ...
    'sltoolbox:slpoolcov:invalidarg', ...
    'w should be a numeric matrix with n columns.');

K = size(w, 1);

if nargin < 3
    vmean = [];
elseif ~isempty(vmean)
    assert(isnumeric(vmean) && ndims(vmean) == 2 && isequal(size(vmean), [d, K]), ...
        'sltoolbox:slpoolcov:invalidarg', ...
        'vmean should be a numeric matrix of size d x K');
end

%% main

% compute mean

if isempty(vmean)
    vmean = slmean(X, w);
end

% normalize w

w = w * (1 / sum(w(:)));

% compute variances

txx = sum(bsxfun(@times, sum(w, 1), X .* X), 2);
tuu = sum(bsxfun(@times, sum(w, 2)', vmean .* vmean), 2);
txu = sum((X * w') .* vmean, 2);

V = txx + tuu - 2 * txu;



